﻿// Learn more about F# at http://fsharp.net

module Gaussian

open System
open System.IO
open System.Collections.Generic
open MathNet.Numerics
open MathNet.Numerics.FSharp
open System.Diagnostics
open MathNet.Numerics.LinearAlgebra.Double

//set point as default decimal separator
System.Threading.Thread.CurrentThread.CurrentCulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")
System.Threading.Thread.CurrentThread.CurrentUICulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")

type IrisData =
  {
    X : DenseVector
    Y : float
  }

let test = File.ReadAllText("iris.data")
let dimensions = 4 //dimensions of the X vector
let dataList = test.Split(',','\n') |> Seq.toList |> List.filter(fun s -> s <> "")
let iris_vector = ResizeArray()
let seed = Random()

let rec gather_data data counter x_vector =
  if data = [] then
    ()
  else
    if counter = dimensions then
      let data_record =
        {
          X = DenseVector(x_vector |> List.toArray) //vector containing the training set input data
          Y = match data.Head with
              | "Iris-setosa" -> 0.0
              | "Iris-versicolor" -> 1.0
              | "Iris-virginica" -> 2.0
              |_ -> -1.0
              //vector containing the expected output
        }
      iris_vector.Add(data_record)
      gather_data data.Tail 0 []
    else
      gather_data data.Tail (counter + 1) (x_vector @ [data.Head |> Double.Parse])

gather_data dataList 0 []

Console.WriteLine("Enter the number of elements used as training set chosen randomly among the Iris set data (negative value means
all elements")
let training_size = Console.ReadLine() |> Convert.ToInt32
let whole_set = training_size <= 0 || training_size >= 150

if training_size <= 0 || training_size >= 150 then
  Console.WriteLine("You either chose an amount less or euqal to zero or greater than the data size. Using whole data set as training
  set. Press a key to continue")
  Console.ReadLine() |> ignore


let pick_random_subset (list:List<'a>) (res:List<'a>) size =
  for i = 1 to size do
    let picked_element = list.[seed.Next(0, list.Count)]
    list.Remove(picked_element) |> ignore
    res.Add(picked_element)
  res
    
let training_set = 
  if whole_set then 
    iris_vector
  else
    let random_set = ResizeArray()
    for i = 0 to training_size - 1 do
      random_set.Add(iris_vector.[i])
    pick_random_subset random_set (ResizeArray()) training_size

let indicator_function element value =
  if element = value then 1.0 else 0.0

let compute_phi class_value =
  (training_set |> Seq.sumBy(fun data -> indicator_function data.Y class_value)) / (float training_set.Count)

let compute_u class_value =
  let class_list = training_set.FindAll(fun data -> data.Y = class_value)
  let mutable res = DenseVector([|0.0 ; 0.0 ; 0.0 ; 0.0|])
  for data in class_list do
    res <- res + data.X
  res <- res / (float class_list.Count)
  res

let phi0 = compute_phi 0.0
let phi1 = compute_phi 1.0
let phi2 = compute_phi 2.0

//DEBUG
//Console.WriteLine("PHI0: " + phi0.ToString())
//Console.WriteLine("PHI1: " + phi1.ToString())
//Console.WriteLine("PHI2: " + phi2.ToString())

let u0 = compute_u 0.0
let u1 = compute_u 1.0
let u2 = compute_u 2.0

//DEBUG
//Console.WriteLine("U0: " + u0.ToString())
//Console.WriteLine("U1: " + u1.ToString())
//Console.WriteLine("U2: " + u2.ToString())

let sigma =
  let mutable res = DenseMatrix(training_set.[0].X.Count)
  for data in training_set do
    let class_mean  =
      match data.Y with
      | 0.0 -> u0
      | 1.0 -> u1
      | _   -> u2

    res <- res + (DenseVector.OuterProduct((data.X - class_mean),(data.X - class_mean)))
  res <- res.Divide(float training_set.Count):?>DenseMatrix
  res

//Console.WriteLine("SIGMA: " + sigma.ToString())

let d = sigma.Determinant()
let gaussian (x:DenseVector) (u:DenseVector) =
  let b = 1.0/(Math.Pow(2.0 * Math.PI,float training_set.[0].X.Count / 2.0) * Math.Pow(d,1.0/2.0))
  let e = (-1.0/2.0) * (x - u) * (sigma.Inverse()) * (x - u)
  b * Math.Exp(e)

let mutable miss = 0
let (output:List<IrisData*float>) = ResizeArray()

let build_output() =
  for data in iris_vector do
    let pxy0 = (gaussian data.X u0) * phi0
    let pxy1 = (gaussian data.X u1) * phi1
    let pxy2 = (gaussian data.X u2) * phi2

    //DEBUG
  //  Console.WriteLine("PXY0: " + pxy0.ToString())
  //  Console.WriteLine("PXY1: " + pxy1.ToString())
  //  Console.WriteLine("PXY2: " + pxy2.ToString())

    let m = [pxy0;pxy1;pxy2] |> Seq.max
    let mutable detected_class = -1.0

    match m with
    | x when x = pxy0 -> detected_class <- 0.0
    | x when x = pxy1 -> detected_class <- 1.0
    | _ -> detected_class <- 2.0

    output.Add(data,detected_class)

let print_output() =
  for (data,detected_class) in output do
    
    let flower_string flower_class =
      match flower_class with
      | 0.0 -> "Iris-setosa"
      | 1.0 -> "Iris-versicolor"
      |_ -> "Iris-virginica"

    if data.Y <> detected_class then
      miss <- miss + 1
      Console.WriteLine("Vector [" + data.X.ToString() + "] has been missclassified as " + flower_string detected_class)
    else
      Console.WriteLine("Vector [" + data.X.ToString() + "] has been classified as " + flower_string detected_class)

let timer = new Stopwatch()
timer.Start()
build_output()
timer.Stop()
print_output()
Console.WriteLine("Missclassifications: " + miss.ToString())
Console.WriteLine("Elapsed time: " + timer.ElapsedMilliseconds.ToString() + " ms")